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1. Introduction 

The experimental progress in the field of ultracold atoms in periodic lattices |3] allowed 
for a direct observation of quantum phase transitions. A very prominent example is the 
experiment by Greiner et al. |4] , where the Mott insulator-superfiuid [5] phase transition 
of ultracold Rubidium atoms trapped in a three-dimensional optical lattice has been 
observed. This phase transition is driven by two counteracting factors. On the one 
hand there is a repulsive interaction between particles on the same lattice site, trying 
to minimize the particle number fiuctuation per site and on the other hand there is a 
gain of kinetic energy whenever a particle tunnels from one site to another. 

Consider a condensate of ultracold atoms in a Mott insulator phase, i.e. every 
lattice site is occupied by the same number of particles. One can then drive a phase 
transition to the superfiuid phase, where the wave functions of the particles are no longer 
localized, by gradually decreasing the intensity of the lattice forming laser beams. 

Because of this tunability of the experimental system [6] it provides a great testing 
ground for various quantum many body theories. Two such theoretical models are the 
Bose-Hubbard (BH) and the Jaynes-Cummings lattice (JCL) model. While the BH 
model considers the competing interaction of kinetic and potential energy of bosons 
explained above |^, the JCL model takes this thought further and also factors in the 
interaction of bosons with a two level atomic system. It therefore constitutes a pure 
quantum transcription of the semi-classical Rabi model. 

Strong coupling approaches have already been used to study the BH model in great 
detail O HI [101 [HI [12] • In this work we treated BH and JC systems with a diagrammatic 
approach to the perturbation theory based on the Kato-Bloch expansion as suggested by 
Eckhard et al. [1] and already applied to the ordered BH model by Teichmann et al. [2]. 

We extended the strong coupling Kato-Bloch approach to deal with disordered BH 
systems and the JCL model with additional degrees of freedom per unit cell. In addition, 
we give details for an efficient numerical computation of the strong coupling diagrams 
and also outline the limitations and possible pitfalls of this approach. 
This paper is organized as follows: The first part deals with the BH model. In section [2] 
we are going to introduce the BH model in more mathematical detail. Section [3] then 
goes on to explain how the Kato formalism can be incorporated to calculate corrections 
to the ground state energy. After this introduction we are going to discuss how this 
approach can be turned into numerical algorithms in section 13.21 

Section 13.31 then shows a comparison of results calculated with the strong coupling 
approach with data obtained with the variational cluster approach (VCA) [13]. In 
section H] we are going to discuss the Mott insulator-superfiuid phase transition, how 
the phase boundary can be detected and which changes have to be applied to the 
present algorithms. Sec. 14.2.21 contains the results we obtained and also a discussion of 
the problems one encounters when trying to determine the phase boundary for a one- 
dimensional system within this approach. In section 14.31 we offer a workaround for this 
problem and show a comparison of our data with results obtained with a density matrix 
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renormalization group (DMRG) approach by Kiihner et al. [H]. Section O explains how 
disordered BH systems can be treated with Kato-Bloch approach. 
In the second part of this work we are going to deal with the JCL model. In section [6] 
we will offer a very short introduction to the JC-model and section 16.11 then outlines 
the numerical transcription and what has to be changed in comparison to the BH 
implementation. The obtained results are shown in section 16. 2[ where we compare 
the energies of various systems for different occupation numbers and dimensions. 
Additionally we compare data from the Kato-Bloch approach with results obtained 
with VGA p5]. The calculation of the Mott insulator-superfluid phase transition is 
explained in section 16.31 

We conclude this work in section [7] with a summary of the obtained results, insights and 
a short outlook. 



2. The Bose-Hubbard model 



The Bose-Hubbard model is used to describe bosonic particles in a lattice at very low 
temperatures [16] and has been treated with many different methods, including the 
mean-field method [7], strong-coupling approaches P El E], DMRG [T71 p:8| [T9| [2n| [21] 
and quantum monte carlo (QMG) methods [22l |23| |2U |25]. A particularly efficient 
strong coupling approach, based on the Kato-Bloch perturbation theory, has been used 
by Teichmann et al. [2]. In the following we will present in short all the definitions 
necessary to understand the subsequent equations. For a more detailed description the 
reader is referred to the above mentioned works. 

The BH hamiltonian H, which we split into two terms suitable for strong coupling 
perturbation theory, reads 



H =Ho + Hi (1) 

Ho := —'^ni{ni - I) - ii'^hi -'^SiUi (2) 

i i i 

H^:=tJ2hb- (3) 



Ho describs the atomic limit, while Hi covers the hopping processes. U stands for 
the strength of the local Goulomb interaction, fi represents the chemical potential in the 
grand canonical ensemble, £j gives the strength of the disorder on site i, n is the particle 
number operator and a and are the well known bosonic annihilation and creation 
operators. The indices label the sites of the lattice. The sums in Hq run over all lattice 
sites, while the sum in Hi runs over all directed bonds (dibonds) hh of nearest neighbour 
sites, i.e.: 

^b = «l%. (4) 
with beV:= {(Z2, ^i)Ki/2 G {1, 2, . . . , iV}} . 

To begin with, we consider a homogeneous system = 0. Disorder can easily be 
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accounted for in a second step as shown in section O Energies will be expressed in units 
of U. The eigenvectors of Hq are occupation number tensor products 

N 

|n) := \ni)i , 

i=l 

where stands for the total number of sites and the number of bosons at site i is given 
by the integer n^. The corresponding eigenvalues are 

N ^ 

■=^n,(^-{ni-l) - (5) 

i=l 

The ground state of Ho is denoted by |g) = <^i\g)i, where each site is occupied by the 
same number of bosons, g, which is determined by minimizing e(n) for a given value of 
the chemical potential, leading to a ground state energy eo = (g|-ffo|g)- 

3. Ground State Energy 

In order to calculate the ground state energy in strong coupling perturbation theory, we 
make use of the Kato-Bloch formalism [31 [2S1 [271 EB] , which yields closed expressions 
for every order of the perturbation, in contrast to standard Schrodinger-Rayleigh 
perturbation theory. The strong coupling Kato-Bloch approach is described in great 
detail in refs. [HE]. 

Here we summarize the key points of the formalism necessary to understand the 
generalization towards disordered systems and the JCL hamiltonian. In addition we 
will present a different perspective that allows to exploit graph theoretical techniques in 
order to speed up the algorithmic implementation. In the form presented in the current 
section, the Kato-Bloch perturbation theory is only applicable to a non-degenerate 
ground state, but we will lift that constraint in section 14.31 
The n^^ order correction for the ground state can be written as 

A4")= Y1 {s\HiSk„_,...H^Sk,H,Sk,H,\g), (6) 

where the sum runs over all sets {fcy}*^""^^ := {(/ci, /c2, • • • , /^n-i)} of non-negative 
integers k^, which satisfy the following conditions [28] 

s 

ki < s ; for s = 1, . . . , n — 2 (7) 

1=1 

n-l 

=n-l. (8) 

1=1 

The operators Sk are diagonal in the occupation-number basis 

(n| Sk |n') = 6n,n' Sk{n) (9) 
-(5n,g for /c = 

Skin) = { 1 - ^n,g . (10) 

r otherwise 

(eo - e{n) f 



Strong coupling expansion for the BH and JC lattice model 



5 



Eventually, the total energy correction can be expressed as a series in powers of the 
perturbation strength t with expansion coefficients a^" 

n 

3.1. Graphical representation 

Next we will map the evaluation of ([6]) onto a graphical problem. For simplicity we 
consider simple cubic lattices in d dimensions. The generalization to other lattices is 
fairly straight forward. We start out from a finite cubic lattice with N lattice sites 
enumerated in some suitable way and periodic boundary conditions. The energy per 
site is computed for the thermodynamic limit, N ^ oo. We consider (|6]) for a specific 
order n and sequence {k,y}^"'~^\ From the n factors Hi we obtain a multiple sum over 
dibonds 

= r J2 (gi ig) • (12) 

There are two important aspects worth noting. Firstly, the application of any of 
the dibond hopping operators hb on an occupation- number basis vector results in just 
another such basis vector 

hb |n) = |n') . 

Secondly, as the operators are diagonal in the occupation number basis they merely 
introduce weight factors. Hence, the index r of the dibond operators hbM defines an 
auxiliary time and 6^"^^ can be any of the 2d dibonds on the d-dimensional simple cubic 
lattice. Each of the individual dibond hopping processes can be depicted as a directed 
line (vector) on the underlying lattice, connecting neighbouring sites (see for example 
figure [T]). The resulting pattern of arrows can be considered as a labelled digraph [29]. 

So far, the sum over dibonds in ( IT2l) can be replaced by a sum over all labeled 
digraphs Q, consisting of directed nearest neighbour lines. Such a digraph covers a 
certain set of lattice sites, which in ordered form be 

ii < i2 < . . . < im , 

where m is the number of vertices of the digraph. Due to the homogeneity of the lattice 
we can renumber the lattice sites in sequential order without changing the resulting 
matrix elements in ( IT2l) . 

iy u ] z/ = 1, . . . , m . 

The digraph defines a set S = {6^*-*} of dibonds, but not the order in which they occur 
in ( !T2|) . Therefore, in addition to the sum over all digraphs Q with n directed lines, we 
have to sum over all distinct sequences of dibond operators occurring in Q. Two such 
distinct sequences of the same digraph are depicted in figure El Obviously, multiple 
dibonds are indistinguishable and need not be permuted. 
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For each dibond sequence b, we compute the corresponding matrix elements 
(weights of the digraph) 

= (g| ht,(n)Sk„_, . . . hh(2)Sk,hhW |g) 

= {gi, ■ ■ ■ ,gm\KwSk„_^ ■ ■ ■h^c2)Sk^h^m\gi, ■ ■ ■ ,gm) , (13) 

where, due to the tensor structure of the vectors and the structure of the operators, 
only the sites 1, . . . , m are involved. An immediate consequence is that the weight of all 
isomorphic digraphs is the same. Since the computation of the matrix elements is the 
most cpu-time consuming step, it is expedient to represent the digraphs by topologically 
different members and the corresponding multiplicities. 

3.1.1. Translational invariance One class of isomorphism is due to translational 
invariance. Each digraph, apart from the lattice-site labels, occurs in N copies on 
the lattice with periodic boundary conditions, which results in an overall factor A^. In 
order to avoid the computation of N identical contributions, we compute the energy per 
site, by restricting the sum over all digraphs to those attached to the origin. Graphs 
attached to the origin shall be defined by the condition 

minf(^)=0, (14) 

i 

where a;'-*^ is the position of the i-th. site on the underlying lattice, covered by the digraph. 

3.1.2. Additional constraints on the graphs Yet not all digraphs contribute to the 
energy. There are two constraints for non-zero contributions. Firstly, the number of 
dibonds entering a site has to balance the number of dibonds leaving it {local particle 
conservation rule). In graph-theory language the digraph is said to be balanced and 
consequently it has no sources or sinks, i.e. sites with only one dibond attached to it. 

A second constraint is due to the linked cluster theorem. As shown in ref. P, [2] and 
references therein, contributions of disconnected graphs cancel. A simple argument is 
given by the extensivity of the energy which has to be proportional to A^. Disconnected 
graphs yield higher orders in A^, since each disconnected part can be placed roughly 
speaking on A^ different sites. 

In summary, admissible digraphs are connected and balanced. From graph theory 
it is known [29j that such graphs are Eulerian, i.e. they can be traversed along the 
directed lines by visiting each directed line exactly once. This very feature will be used 
to generate all digraphs. Examples of several low-order diagrams are given in refs. 

3.2. Algorithm 

Here we give an algorithm to generate the topologically different, connected and 
balanced digraphs on a simple cubic lattice. 
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Figure 1. Two identical digraphs of order n = 12. The generating paths are attached 
to the origin. 



3.2.1. Generation of topologically different digraphs We start at the origin of the 
lattice and generate an Eulerian circuit of length n by tracing a tour sequentially to 
neighbouring sites. At every site there are 2d choices to proceed. The problem is 
equivalent to searching in trees and we applied a depth-first search algorithm up to 
depth n. If we return to the origin after these n steps, the corresponding path is shifted 
such that it fulfills and is added to a list. Multiple entries of the same path are 
discarded. The same path can occur several times, as we can consider any site of the 
path as origin, from which the path is initiated. This multiplicity is an artifact of the 
algorithm and has to be eliminated. 

Next we need to identify topologically different digraphs. Digraphs are isomorphic 
if and only if for some ordering of their vertices, their adjacency matrices [30] are equal. 

Let {^l,^2, • • • ,^n+i} be the sequence of site-indices encountered during the tour. 
The adjacency matrix, which includes the information which sites of a graph are 
connected with each other, can be computed in the following way 

n 
1=1 

where Mij is the number of times the path goes from site i to site j or in other words, 
the number of dibonds from site i to site j. 

Two digraphs are isomorphic if and only if the corresponding matrices are similar 

M' = PMP-^ . 

with P being a permutation matrix pO]. Provided two graphs are isomorphic, trace and 
determinant of M'^ and M'" for any integer power u have to agree. This can be used as 
a quick tests to rule out isomorphism. 

A representative example of two digraphs that pass all these tests is depicted 
in figure [H The corresponding adjacency matrix of the left digraph is 
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(15) 



(16) 



At first glance, the two matrices seem to be different, but the distinction stems from 
the labehng. If we denote by P the matrix associated with the permutation 

((1), (2, 4), (3), (5), (6, 7), (8, 9)) , 

then we can easily verify the similarity of M\ and 

M2 = p-^MiP . 

This algorithm then generates a hst of unlabeled (topologically different) digraphs, along 
with their multiplicities. 



3.2.2. Computation of the energy correction In order to evaluate the energy correction 
per site of (|6]) we have to go through all digraphs in the list, compute their contribution to 
the matrix element and multiply it by the multiplicity. The contribution of one digraph 
is determined in the following way. The digraph Q defines a set of dibond operators. 
We attach auxiliary times r = 1,...,?t, to the dibond operators in all possible but 
distinguishable ways. Each time sequence {h^'^\ . . . , h^^^) defines an additive contribution 
to ©. 

There is an implicit condition due to the fact that 

(g| |g) = (17) 

for every dibond operator. Consequently, /ci 7^ and kn-i 7^ and for the same reason, 
consecutive /c-values cannot be zero simultaneously, i.e. ki /cj+i 7^ 0. 
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Figure 2. Two possible time resolved dibond sequences of a digraph of 4'^ order. 



3.2.3. Summary of the algorithm 

• Generation of paths: First we generate all closed paths of a given order n, beginning 
at the origin. 

• Translational invariance: Each path is attached to the origin as defined by fn3|) . 
Multiple occurrences of the same path are eliminated. 

• Digraphs: The paths are translated into digraphs. Site labels are omitted and 
(topologically) unique digraphs are stored along with their multiplicities. 

• Dibond sequences: For each digraph, all distinct time dependent dibond sequences 
are generated. 

• Kato-Bloch summation: For each dibond sequence, the compatible Kato-Bloch 
indices {kv}^'^~^^ are determined and the contribution to the energy correction ([6]) 
are computed. Care is taken when admissible subgraphs occur. 



3.3. Test result 

In order to test the numerical algorithm we compared the ground state energy as a 
function of the hopping strength for the one dimensional BH model with the results 
obtained by VGA [13]. An example is given in figure [3] and as one can see, the results 
are in very good agreement. More results for the BH model obtained with the strong 
coupling Kato-Bloch approach can be found in ref. [2]. 



4. Mott insulator-superfiuid phase transition 

We will now address the phase boundary between the Mott insulator and the superfiuid 
phase, which can easily be computed in the frame of the strong coupling Kato-Bloch 
formalism. There are essentially two approaches, for fixed ii or for fixed t. 

The latter is driven by charge fiuctuations and the phase boundary can be detected 
by computing the energy Ai?^^ it takes to add or subtract a 'defect' particle as has 
been shown in ref. [8] (see section l473l) . 

At fixed /i and density, the transition is driven by phase fiuctuations. To detect 
this transition, we introduce an 'Effective Potential' [H [HI |32] • 
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Figure 3. Energy corrections with the present approach (up to 8*'' order) for 1 
dimension and filling factors 5 = 1 and g = 2 compared with results obtained with 
VGA [13]. 



4-.1. Phase Boundary Criterion 

As a very detailed description of the method of the 'Effective Potential' is given in 
ref. [2], we will do without a derivation and just present the necessary equations and 
definitions. 

In order to describe the superfluid phase, we have to add a source and a drain term to 
the BH hamiltonian 

H = HoJ^Hf^+Y,{v*ai + val) , (18) 

where t] is the strength enforcing particle number fluctuations. The order parameter 
= {0'i)'ri can then be written as a power series in rj 

ij = {c2 + 2ci\7]\^ + 0{\r]\^))r] . (19) 

The coefficients C2n that appear in (fT9|) are expanded in power series of the hopping 
parameter t: 

oo 

C2, = 5^aJ)r (20) 

n=0 

As pointed out in ref. [2], C2 can be considered as a susceptibility 

drj J ri^o ' 

of the system to develop an order parameter driven by the source-and-drain term of 
strength r]. Hence, the Mott-to-superfiuid phase transition occurs at those points in the 
H — t— phase diagram, where C2 diverges. 



C2 
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The convergence radius of the power series of the function Og^^ (t) in ( l20|) is according 
to d'Alembert's law [33] given by 

(21) 



t* = hm 



Or, 



Hence, for given chemical potential, t* marks the transition from the Mott insulator to 
the superfluid phase. 

In order to use the strong coupling Kato-Bloch formalism the perturbation term 
has to be modified 



In analogy to ([HD and f[T^ we have 

«r'^= E E (g|5("^5,_...6(^^5.,6«|g) • (22) 

{fc^}("-i) 6(1), 

The operators 6*^'^^ are either dibond operators hf,(T) as in ([3]) or individual particle 
creation a| or annihilation operators a, at site i. As a matter of fact, since we are only 
interested in the term C2, exactly one creation and one annihilation operator have to be 
present. 

4-2. Graphical representation 

For the n*^ order contribution [a^'^) the graphical elements are n nearest neighbour 
directed lines, one source term and one drain term. By virtue of the linked cluster 
theorem the graph has to be connected and the source and drain symbols are located 
on the vertices of the graph. On top of that, particle conservation demands 

d+(z) -rf"(z) = iV+(z) -iV^(z) , 

where d~^^~{i) are out-/in-degrees of vertex i and N'^^~{i) is the number of source- 
/drain-symbols of vertex i. In other words, vertices with no source/drain or with a 
source-drain pair are balanced, while a vertex with only a source/drain symbol has to 
have one remaining, not compensated outgoing/incoming directed line. 

4-2.1. Construction of admissible digraphs The admissible digraphs can easily be 
completed to form a balanced connected digraph, by adding an extra bath site, from 
which a directed line points to the drain symbol and a directed line pointing from the 
source symbol. This is not necessarily a nearest neighbour line and can even be a self- 
loop. The sought-for digraphs are therefore Eulerian and can again be constructed by 
tracing a continuous tour. The tour begins at the origin with a source (•) and proceeds 
successively along nearest neighbour dibonds. For the n*^ order term ag"^ '^g generate n 
directed lines and close the graph with a drain (x). (See for example figure H] or ref. [2].) 

In contrast to the energy calculation however, it is not compulsory that the path is 
closed. As before, the digraph is attached to the origin according to (fT4|) and multiple 
copies of the same digraph are discarded. 
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Figure 4. Four examples for paths with one source (•) and one dram (x). 

Table 1. Number of topologicaUy different digraphs for given order and dimension 



order 


1-dim 


2-dim 


3-dim 


1 


1 


1 


1 


2 


2 


2 


2 


3 


4 


4 


4 


4 


8 


10 


10 


5 


14 


22 


22 


6 


25 


58 


58 


7 


45 


140 


140 


8 


79 


390 


394 



The expectation value in ( ]22|) is invariant against relabeling the sites, i.e. isomorphic 
digraphs yield the same contribution to a'^\ We again identify isomorphism by similar 
adjacency matrices. 

In addition, table [T]shows the number of topologicaUy unique diagrams as a function 
of the order for systems with 1, 2 and 3 dimensions. The reason, why the number of 
topologicaUy unique diagrams is the same for all dimensions up to third order is that 
with a maximum number of 3 bonds one can only draw paths that can be mapped to 
1-dimensional ones (see (a), (b) and (c) in figure Hj). Only in fourth order or higher is 
it possible to draw a closed loop in two and three dimensions, which is not possible in 
one dimension (figure SU^d)). The same reasoning applies to the 2- and 3-dimensional 
case, which have the same amount of topologicaUy unique paths up to seventh order. 
Only in eighth order or higher paths in 3 dimensions exist that don't have topological 
equivalents in 2 dimensions. 

As before, we have to sum over all distinguishable 'time'-sequences of the elements 
in the digraphs, which are this time n dibond hopping terms, one creation and one 
annihilation operator. 

Following (12T]) . we calculate the ratio |«2'^~^V'^2'^^I a function of l/u. The 
extrapolation ^ — > yields the sought-for value of the phase boundary. In figure \5\ 
a representative example is depicted. Further details and examples can be found in 
refs. [HE]. 





V 

X 
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Figure 5. The ratios \a"2 /ct'^^'l as function of 1/;^ and extrapolated to — >■ oo 
using a linear fit for the first Mott lobe {g = 1) for a 2d system with n/U = 0.3. 



4-2.2. Test results The phase diagram for the BH model has already been computed 
with high accuracy by various techniques as well as by strong coupling approaches [21 
EH El |9] . Our results agree very well with those results. While the 2d and 3d results 
agree remarkably well with results from other methods |35], the Id case causes problems 
close to the tip of the Mott lobe. In figure [6] we compare our results with data from 
DMRG calculations [H]. 

Teichmann et al. [2] already pointed out that the strong coupling Kato-Bloch 
approach has some problems in detecting the phase boundary in Id. They put the 
blame on the reentrance feature of the Mott lobe: For a given chemical potential, the 
criterion based on the radius of convergence in ( l2Ti) can only provide the phase boundary 
with the smallest t*. We observe however, that deviations also occur for values of the 
chemical potential, where there is no reentrance feature. The agreement with the DMRG 
results is very good for small values of t, but close to the tip deviations start to grow. 

The discrepancy in the vicinity of the tip is due to the power-law decay of the 
correlation functions of the Kosterlitz-Thouless (KT) transition at the tip. 

The treatment of such a transition is very difficult as much depends on the 
extrapolation scheme. In a very recent paper, Ejima et al. [21] calculated the tips 
of the first two Mott lobes using DMRG and a suitable treatment of the correlation 
functions. The results are remarkably close to those values obtained with QMC. 

4.3. Method for 1- dimensional systems 

As pointed out in the previous section, the power-law behaviour of the KT transistion 
causes problems for the method of the effective potential. We therefore incorporated 
another approach that allows us to determine the phase boundary for Id systems, 
proposed by Freericks et al. [8]. 

The idea is to introduce defect states, i.e. one adds/removes one particle from 
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the Mott ground state |g) — )■ |. . . ,gi^i,gi ± . . .). With the assumption, that 

the compressibihty approaches zero countinously at the phase boundary, the critical 
chemical potential fi* can be deduced at the point where 

AEp^ = AE^^^ - AE^ = . (23) 

This puts us in the advantageous situation that we can now use the energy algorithm 
described in section 13.21 which has the benefit of being immensly faster: In contrast to 
the method of the effective potential, we now have to consider only closed paths. The 
computation time can therefore be reduced by a considerable amount. Additionally, 
instead of having to calculate the critical t* for every /i of interest seperately, we only 
have to compute AE once for every particle (hole) band and can deduce from this all 
/i* with a simple numerical operation. 

The precise procedure works as follows. We start with the already well known Mott 
insulator ground state |g). Then we add (remove) one particle at the site of the origin 
of the considered path (We shall from now on refer to this new ground state as |p(h))). 
Now we let this particle (hole) hop through the lattice according to the diagram and 
all unique permutations of it. One has to pay attention to the fact that the system for 
a given diagram with N involved sites is A^-times degenerate [28]. In other words, the 
state \g ± 1, g, g) leads to the same unperturbed energy eo as the states \g, g g) and 

While in the nondegenerated case, sequences containing factors of the form SohhSo 
could have been omitted (see. ( ITTl) ). this is no longer allowed in the degenerate case. 
Consequently, all Kato sequences that include a factor (p(h)| hb |p(h)) can indeed have a 
nonvanishing contribution and have to be dealt with as well. Apart from that, everything 
works as explained in section [31 



Strong coupling expansion for the BH and JC lattice model 



15 




Figure [7] shows a comparison of results calculated up to 9*^ order with data from 
Kiihner et al. [H]. As one can see, the perturbation approach results match those from 
DMRG quite well, but deviate discernably for larger t, especially in the particle band 
case (i.e. the upper branch in the phase diagram). 

This is due to the fact that the power series expansion of /i* has an asymptotic 
behaviour. Calculating higher orders of the perturbation series in order to improve 
accuracy is therefore not feasible. Hence we are going another route and do the following: 
We start with the power series for /x* for the particle case (the hole case can be treated 
analogously), which we get from the Kato-Bloch energy calculations 



t^lit) = 2^ , (24) 

and perform a Borel transformation (see for example ref. |43] ) 

B{z) = Y^c,{vf-,- (25) 

This Borel transform is now rewritten and approximated according to Fade in order 
to cope with the asymptotic series. 

^(-) = • (26) 

1+ ^nX" 
n=l 

For a 1-dimensional system with filling factor g = I for example we get the following 
series Bp{x) for the particle case and Bh{x) for the hole case (All coefficients have been 
rounded to two decimal figures): 

1 - 3.11X - 3.21x2 + 1.88x3 + 1.26x^ 
~ 1 + 0.89x - 0.14x2 - 0.15x3 + O.OQx^ 
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Figure 8. Comparison of our results up to 9* order obtained with the Kato approach 
and DMRG from Kiihner et aL [14] for a 1-dimensional system with filhng factor g = 1 
after the Borel-Pade transformation has been performed. 



2x- 1.75x2 - 0.02x3 - 1.27x^ 
^^^^ ~ 1 + 0.13x + 0.12x2 - 0.10x3 + 0.04x4 ' 
The final series for fi*(t) is then recovered by evaluating the integral 

oo 

m-jsMK'^. (27) 



This revised expression for /i* is very good natured and our results now agree excellently 
with those from the literature [H] , as can be seen in figure [H 



5. Local Disorder 



In a previous work, Freericks et al. used Schrodinger-Reighleigh perturbation 
expansions up to third order to study a disordered BH system and Krutitsky et al. [12] 
further investigated such systems with high order strong coupling expansions. We will 
demonstrate how local disorder, described by the last term in ([2]) (where Ei is a random 
variable), can be dealt with in the framework of the cluster Kato-Bloch formalism. 

The effectiveness of the Kato-Bloch approach relies on reducing all paths to just 
topologically unique ones, which is only possible for a translational invariant system. 
By introducing a random disorder the invariance for a single path is broken, therefore 
the central quantity now has to be the grand canonical potential at zero temperature, 
or rather the lowest energy in Fock space for a given chemical potential, averaged over 
disorder realizations. By averaging over a large number of disorder realizations, the 
translational invariance is reobtained. 

{E)^ := j E{e) p{e) d^'e , (28) 

where E{e) is the lowest grand canonical energy for a given disorder realization vector e 
and p{e) is the joint probability density function (pdf ) for the random disorder variables. 
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The latter are assumed to be uncorrelated and the joint pdf is the product of independent 
and identically distributed (iid) random variables 

N 

1=1 

The n}^ order contribution to the mean energy is according to (Ej) 

(^4"^>.= E {is'\HiSl-.---SlH^\gn)^ , (29) 

where the ground state occupation g"^ depends on the disorder realization and will be 
site dependent. The same holds for the energies in the atomic limit E^{n), entering the 
operators Sf.. The summation over all dibond sequences in (fT2l) is again mapped into a 
sum over all connected and balanced labelled digraphs. Similar to ( 1T3|) . the weight for 
one digraph Q is given by 

w{g) = (30) 

where only the disorder energies of the site reached by the graph are required. The 
mean weight is again independent of the labeling and hence isomorphic digraphs yield 
the same contribution. So the algorithm is almost the same as before in the homogeneous 
system. We generate all topologically different, connected and balanced digraphs along 
with their multiplicities. For each digraph all distinct permutations of the dibonds are 
generated and the summation over possible Kato-Bloch sequences is performed. 

For every specific dibond sequence and every set of Kato-Bloch indices the averages 
over the disorder realizations on the m involved sites has now to be carried out. The 
same modification holds for the computation of the susceptibility a^- 

In order to compare with ref. [36], we have chosen a binary disorder with Si = ±e, 
with disorder strength e. 

5.1. Results 

Figure |9] shows the energy for a 1-dimensional system as function of the hopping 
parameter for different disorder parameters e as well as a comparison with results 
obtained with VGA [36]. Energy corrections were included up to 8**^ order and we 
averaged over 500 configurations. The solid lines are our results calculated with the 
Kato-Bloch algorithm and the x mark the results from VGA. 

As before, our results are in excellent agreement with those computed with VGA, 
which can be seen especially well for the e = 0.1 and e = 0.2 cases. Unfortunately we 
don't have data from VGA for e = 0.3 and e = 0.4 at higher hopping parameters, so 
the possibilities to compare these two methods for these parameters are rather little. 
The data at low hopping strength and the general information gathered up to this point 
however suggest a good agreement of Kato-Bloch and VGA nevertheless. 
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Figure 9. Energies of a 1-dimensional system with disorder parameters e = 0.1 (blue 
dotted line), e = 0.2 (red dashed-dotted line), e = 0.3 (green dashed line) and e = 0.4 
(black solid line) compared to results obtained with VGA [35] , which are marked with 
X. The chemical potential is fixed at /i = 0.5. 




Figure 10. The Mott insulator gap for a 2-dimensional system as function of the 
hopping strength for e = 0.0, e = 0.1 and £ = 0.2. 



In figure [TD] we plotted the width of the Mott insulator regions as function of the 
critical hopping strength t* as a way for comparison. The blue solid line stands for the 
ordered system, the red one for the system with e = 0.1 and the black one for e = 0.2. 
In figure [TT] the dependence on the disorder parameter e is depicted for the g = 1 Mott 
lobe. We observe that the Mott lobe shrinks with increase disorder. A comparison of 
the 2-d and 3-d reveals that the impact of the disorder is almost independent of the 
physical dimension. It is important to note however, that within the present approach 
it is not possible to tell whether the neighbouring phase is a superfiuid or a Bose glass 
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.04 



r/u r/u 

Figure 11. Mott insulator for g = 1 for a 2-dimensional system (left figure) and a 
3-dimcnsional system (right figure) as function of the hopping strength for e ~ 0.0, 
e = 0.1 and e = 0.2. 

phase. 

6. The Jaynes-Cummings lattice model 

Due to experimental progress in controlling quantum optical and atomic systems, new 
ideas have been proposed for new realizations of strongly-correlated many body systems, 
such as ultracold gases of atoms trapped in optical lattices [HI SI [3] or hght-matter 
systems. ESI EH] The latter contains photons, that interact with atoms or atomic- 
like structures. The interaction is strong if the photons are confined within optical 
cavities. SH Il2] The coupling between photons and atoms leads to an effective 
repulsion between photons, which in turn results physical properties like in the Bose- 
Hubbard model, such as the quantum phase transition from a Mott phase, where 
particles are localized on the lattice sites, to a superfluid phase, where particles are 
delocalized on the whole lattice. [37] Yet the physics of the light-matter models is far 
richer because two distinct particles, namely photons and atomic-like excitations, are 
present. 

The following short explanation follows strongly the approach of ref. [16] and all 
quantities are expressed in units of h. 

The hamiltonian of a JC system consists of three parts, an atomic part Ha = e |t) (tl 
which assigns the energy e to the excited atom states, a cavity part He = WcS^a counting 
the bosons in the cavity and appointing them the energy Uc and a part that describes 
the coupling between atom and cavity Hac, which can be written like 

Hac = + ) (31) 

with the Rabi frequency Q, which is a measure for the strength of the coupling of the 
two systems. Here, a± corresponds to the atomic raising and lowering operators, while 
a and are the usual photonic annihilation and creation operators. 

The eigenstates of Ha are \D, the first denoting the ground state, the second 
the excited state. On the other hand, the cavity hamiltonian He has the eigenstates 
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|n), the already known Fock-states from section [2J Therefore, the eigenstates of the 
uncoupled system will be the tensor products \n,\) and |r;,,t). 

If the detuning of the system A = — e is zero or very small compared to cUc, the 
states with the same particle number, i.e. |n, J,) and \n — are degenerate or nearly 
degenerate respectively. The complete energy of a system with n particles is therefore 
saved in a state with n photons and no atomic excitation \n,\) and a state with n — 1 
photons and one atomic excitation |n — l,t)- (The exception being a system with no 
particles, which of course can only be described by |0,i).) The coupling hamiltonian 
Hac only translates between those states with the same particle number. 

In this case, the energy eigenvalue can be written like 

E„,± = ncuc - ^ ± q{.n) , (32) 



with q[n) = \l \ — + n [ — 



2 J \2 

The ±-sign refers to the sign of \n, — ) and |n, +), introduced in ( 1331] and 0341] . If there are 
no particles in the system, there is just one state |0,i) possible. The energy eigenvalue 
in this case is -E„=o = 0. 

The eigenstates corresponding to the energies of fl32|) are 

|n, — ) = cos 6„, |n — 1, t) — sin 6„ |n, I) (33) 
|n, +) = sin 9„ |n — 1, t) + cos 6„ J,) , (34) 

where we used the short notations 



q{n) - 


A 
2 


2q{n) 


\q{n) + 


A 
2 



sin0„='i/ — ^ and (35) 



cos 9„ 



2q{n) 

Additionally to what was said in the previous paragraphs, an atom-cavity site may also 
have an additional energy which is dependend on the total particle number, due to a 
chemical potential /i. The final JC hamiltonian for an atom-cavity system looks as 
follows: 

Hjc = ujca^a + e\^){^\+ (36) 

+ §(«lt) (t|)-/i(a^«+lt) (tl) • 

Now we are going to construct a regular lattice by arranging a large array of such sites 
and allow the bosons to tunnel to neighbouring sites with tunneling strength t. The JC 
hamiltonian of a single site Hjc,i is given in ( |371) and the hopping hamiltonian is the 
same as for the BH model (see section [2]) 

HjcL = ^ Hjc,i + Hhop ■ (37) 
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Figure 12. One possible path in second order. 
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Figure 13. All possible sequences of the path in figure [ 



6.1. Changes in the algorithms 

The JC site hamiltonian Hjc does not alter the particle number of the eigenstates, but 
the hopping hamiltonian Hhop does. It is therefore clear that the hopping term is again 
the perturbative part. 

The fact, that there are two eigenstates of the JC model now belonging to the same 
particle number n involves a disadvantage. As explained before, the term 5*^ of our Kato 
formula projects the current state onto the according eigenstate of Hq, i.e. it will either 
be projected onto \n, — ) or \n, +). This means, that one has to add the information to 
which eigenstate we are projecting, i.e. 5*^ — )■ Sk,a with a indicating which eigenstate 
we should take: 

— Sn,gSa~ ioY k = 

' otherwise. 

(38) 



S, 



k,± 



^{0) _ ^(0) 



In order to calculate the energy correctly we also have to consider not only every possible 
sequence {ki,Y'^~^^ but also every possible permutation of the signs {cry}^""^^. 
Our equation for the complete energy correction of n*^ order therefore reads 

AeW= J2 ^^il- (39) 

To deduce the right sequences we will take a look at a few different paths. The most 
simple path appearing in the second order energy correction is depicted in figure [121 
After the first hop from site 1 to site 2, we do have g — 1 particles at site 1. But the state 
of this site could have either been projected to \g — 1,—) or \g — 1,+). An exception 
would of course be if there were g = 1 particles at each site before the hopping, in which 
case only the \0,l) state would be possible for site 1 after the first hopping. We will 
from now on omit mentioning this exception as it always results in the same state. 
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Figure 14. One possible path in fourth order with hopping sequence {1,2-2,3-3,2-2,1}. 



Site 2 contains g + 1 particles after the first tunneling and will be projected either 
onto the \g + 1, — ) or the \g + 1, -|-) state. After the second hopping takes place, both 
sites are again populated by g particles and at each site the ground state \g, — ) has to be 
present in order to lead to a non-vanishing energy contribution. A graphical depiction 
of this can be seen in figure [131 

To calculate the energy contribution of this path correctly, we would now have to 
calculate the energy of each of the 4 possible sequences and sum them up. 

In order to get the right number of sequences, we have to include information about 
the number of nodal points, i.e. how many times the ground state is restored before the 
whole perturbation process is finished. An instructive example is shown in figure [Ml and 
figure Uni While both of these figures show the same diagram, the order of hoppings is 
different. In figure [Til the ground state is never present up until the very end resulting 
in the sequences of figure [151 For the path depicted in figure [161 however we assume 
it restores the ground state after two hopping processes, leading to the sequences of 
figure [T71 (If the path of figure [IHl does not result in the ground state after the first two 
hoppings but in the |g, -|-) state, we would get a very similar diagram for the sequences 
as that of figure [T5l ) To compute the total numbers of sequences for a certain path, 
we count the number of arrows pointing to a site Na,i- For closed paths this number has 
to be of course the same as the number of arrows pointing away from a site, the number 
of bonds of site i is therefore Ns^i = 2Na^i. Additionally, we need to know the number 
of nodal points N^^i for this site. The number of sequences Ngeq^i for site i is then 

Nseq, = 2^«--(^"- + l) , (40) 

and the total number of sequences Ns is in further consequence the product of these 



Ns = l[N,eq,^■ (41) 



The number of bonds for a specific site can be determined very easily with the aid of 
adjacency matrices. The i^^ row contains the number of bonds originating from site i. 
The i^^ column on the other hand contains the information from which sites bonds point 
to site i. The adjacency matrix for the path shown in figure is 

/ 1 \ 

M= 10 1 . (42) 

{o 1 oj 

The vector vb resulting when summing over all columns of the matrix stated by (H2l) 
reads = (1,2,1). Doubling it results in (2,4,2), which means that there are two 
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Figure 15. All possible sequences of the path in figure [T4l 



Figure 16. One possible path in fourth order with hopping sequence {1,2-2,1-2,3-3,2}. 



bonds at site 1, four bonds at site 2 and again two bonds at site 3. 

With all that in mind we can reuse most of the algorithms created for the energy 
calculations of the BH model. The path creation and the reduction to topologically 
unique diagrams can be copied without any change. The algorithm for the energy 
calculation however has to be adopted at the point where all the permutations of a 
diagram are created. At this point we have to compute the adjacency matrix for every 
permutation to determine the number of bonds for each site. With this information 
we can create the allowed cr-sequences for each permutation individually, calculate the 
energy for every sequence and sum them up. After all cx-sequences have been worked 
off, we proceed to the next permutation. 
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Figure 17. All possible sequences of the path in figure [TBI 



Table 2. Computation times in seconds for the JCL energy correction up to order 
v — Q for different filling factors g and dimensions d. (All computations were done 
with Matlab on a computer with a clock frequency of 3 GHz.) 







9 = 1 






9 = 2 






d = 1 


d = 2 


d = 3 


d=l 


d = 2 


d = 3 


2 


0,25 


0,24 


0,35 


0,23 


0,22 


0,32 


4 


0,37 


0,46 


0,53 


0,58 


0,72 


0,75 


6 


8,30 


27, 15 


27,71 


58,61 


153,06 


156,54 



6.2. Results 

The additional summation over all a-sequences of course slows down the energy 
computation compared to the calculations for the BH model, but the computation times 
up to order 6 are still in the range of seconds, as table [2] shows. 

All systems we are going to discuss have in common that we set the detuning A = 0. 
A more detailed look at different sets of parameters can be found in ref . [IE] . 

Figure [18] shows the computed energy correction of a 2-dimensional system with 
filling factor g = 1 for different orders of correction. The solid (blue) line contains only 
the second order energy correction AE^'^\ the dashed (red) line results from combining 
the second order and fourth order corrections, the dashed-dotted (green) line shows the 
energy correction up to sixth order and the dotted (black) line up to the eighth order. 

While the solid (blue) line differs from the others quite significantly, the dashed 
and dotted lines are very close to each other. In fact, the lines for the 6*^ and 8*'* order 
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Figure 18. Comparing the fidelity of the computed energy corrections depending on 
the orders of correction included for a JCL system with 2 dimensions and unity filling 
factor. 
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Figure 19. Comparison of calculations carried out with the Kato-Bloch algorithm 
(solid and dashed lines) and with VGA [TS] (circles and x) of two systems with 
dimension d = 2 but different fillings factors g. 



lie on top of each other perfectly at this scale, suggesting a very high fidelity already 
for the sixth order correction. 

Additionally, figure [19] shows a comparison of data computed with VGA [15] . Two 
systems were simulated, with the first having 2 dimensions and unity filling factor and 
the second also 2 dimensions but filling factor g = 2. 

The solid and the dashed lines represent the results obtained with the Kato-Bloch 
approach and the circles and x those from VGA. As previously with the BH system 
(see figure [3] in section [373]) . our results are in excellent agreement with VGA. 
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A comparison of the energy of different systems witli 1, 2 and 3 dimensions and 
with the fiUing factors g = 1 and g = 2 is shown in figure [201 Based on the numerical 
resuhs, it appears that the energies depend stronger on the fihing factor g than on the 
physical dimension d. 

6.3. The Mott Insulator- Superfluid Phase Transition 

The same changes in the algorithm that have been discussed in the previous chapter have 
also to be taken into account if one wants to determine the Mott insulator-superfluid 
phase boundary for the JCL model. But furthermore it is no longer necessarily true 
that the number of bonds pointing to a site is the same as the number of bonds pointing 
away from it. This little inconvenience is taken care of by just computing the sum along 
the rows and the sum along the columns of the adjacency matrix and adding these 
vectors point- wise, thereby getting the number of bonds per site. In Figure [21] the Mott 
insulator-superfluid phase boundary is shown for d = 2 and d = 3 and different filling 
factors. Left to the depicted lines, the system is in the Mott phase. 

The phase boundary Hq for t* = is independent of the dimension of the system as it 
is determined by a single site. It can be calculated easily, resulting in /Xq = ^/g — ^ g + 1 
(see ref. [^). This A comparison with the result for the BH model, where fi^ = Ug, 
reveals that the effective Hubbard interaction U in the Jaynes-Cummings model depends 
on the filling fractions. This explains, why the Mott lobes in the JC model decrease 
rapidly with filling fraction. 

As was the case with the BH model, the phase boundary for 1-dimensional systems 
has to be calculated by the use of 'defect' particles already explained in section [4751 
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Figure 21. Boundary of the Mott phase for d = 2 and d = 3. The detuning is zero 
A = 0. The filling factors for the three Mott lobes are g = (0, 1, 2) from the bottom 
to top. 

7. Conclusion and Outlook 

In this work we have employed the Kato-Bloch perturbation approach to calculate the 
ground state energies and boundaries of the Mott-insulator phase, both for various 
Bose-Hubbard and the Jaynes-Cummings lattice systems. 

The Kato-Bloch approach leads to explicit formulas for any perturbation order, i.e. 
there is a clear distinction of the orders in contrast to the mixed expressions in the 
Schrodinger-Rayleigh perturbation theory. Moreover, all the appearing expressions can 
be represented diagrammatically, which allows to utilize knowledge from graph theory 
to speed up the numerical algorithms considerably. In fact we were able to perform all 
calculations on an off-the-shelf computer using Matlab in appropriate time. 

Additionally to these advantages, the comparison of our results with those 
from other methods, including DMRG and VGA prove the high accuracy of our 
implement at ion . 

Building upon the work by Teichmann et al. [2], who determined the ground state 
energy, the boundary of the Mott-insulator phase and various other quantities for the 
BH system in detail, we focused on adding disorder to the BH model, applying this 
approach to the JCL model and on solving the problem that appears when trying to 
calculate the phase boundary for 1-dimensional systems with the method of effective 
potential. 

The Kato-Bloch approach is an extremely powerful method and could be easily 
generalized to deal with other lattice systems apart from the BH and the JCL as well. 
Especially interesting in our eyes would be the implementation of models with nearest 
neighbour interaction or an additional superlattice structure, which exhibit a richer 
phase diagram as the BH and the JCL model. 
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